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Abstract 

The low-lying eigenvalue spectrum of the QCD Dirac operator in the e-regime is expected to 
match with that of chiral Random Matrix Theory (ChRMT). We study this correspondence for 
the case including sea quarks by performing two-flavor QCD simulations on the lattice. Using the 
overlap fermion formulation, which preserves exact chiral symmetry at finite lattice spacings, we 
push the sea quark mass down to ~ 3 MeV on a 16^ x 32 lattice at a lattice spacing a ~ 0.11 fm. 
We compare the low-lying eigenvalue distributions and find a good agreement with the analytical 
predictions of ChRMT. By matching the lowest-lying eigenvalue we extract the chiral condensate, 
S^^(2 GeV) = (251 it 7 it 11 MeV)^, where errors represent statistical and higher order effects 
in the e expansion. We also calculate the eigenvalue distributions on the lattices with heavier sea 
quarks at two lattice spacings. Although the e expansion is not applied for those sea quarks, we 
find a reasonable agreement of the Dirac operator spectrum with ChRMT. The value of S, after 
extrapolating to the chiral limit, is consistent with the estimate in the e-regime. 
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I. INTRODUCTION 



Numerical simulations of QCD on the lattice suffer from various sources of systematic 
errors, such as finite lattice spacing a, finite volume V, and larger quark masses m than 
those in the nature. Each of these needs to be eliminated by an extrapolation using several 
independent simulations. In particular, the extrapolation in the quark mass to the chiral (or 
physical) limit is non-trivial, because most physical quantities have non-analytic dependence 
on the quark masses due to pion loop effects as predicted by chiral perturbation theory 
(ChPT). In order to reproduce such non-analytic behavior, the physical volume must be 
increased as the chiral limit is approached such that the pion Compton wavelength fits in 
the box. Therefore, in practice the chiral extrapolation must be done with a limited range 
of quark masses, which is a potential source of large systematic uncertainty. This becomes 
more problematic when the chiral symmetry is explicitly violated by the fermion formulation 
on the lattice, since the standard ChPT cannot be used as a guide in the extrapolation and 
the chiral extrapolation must be combined with the continuum extrapolation. 

An alternative approach is to study the e-regime of QCD 

mm 

on the lattice. In this 

regime the quark mass is set close to the chiral limit while keeping the physical volume finite. 
The system suffers from a large finite volume effect, but it can be systematically calculated 
by ChPT, because the pion field dominates the low energy dynamics of the system and 
the effects of other heavier hadrons become sub-dominant. It means that the low energy 
constants appearing in ChPT Lagrangian can be extracted from the lattice calculation in the 
e-regime by comparing with ChPT predictions. Since a small violation of chiral symmetry 
gives large effects in the e-regime, the lattice fermion formulation must fully respect the 
chiral symmetry. 

The e-regime is reached by reducing the quark mass m, at a finite volume V = L^T, 
down to the region where the pion mass satisfies the condition 

1/Aqcd < < (1) 

where Aqcd denotes the QCD scale. Under the condition ([T]), the zero momentum modes 
of the pion field give the dominant contribution since the energy of finite momentum modes 
is too large to excite. In this way, ChPT is organized as an expansion in terms of the 
parameter e^ ~ mT^/Auy ~ p^/A^v where Ayy is the ultraviolet cut-off of ChPT (typically 
taken to be 4:nF^ with the pion decay constant). Since the quantum correction of the 



zero-modes is not suppressed in the e-regime and the path integral over SU(A^/) manifold 
must be explicitly carried out, the partition function and other physical quantities show 
remarkable sensitivity to the topology of the gauge field. 

At the leading order of the e-expansion, the partition function of ChPT is equivalent 
to that of chiral Random Matrix Theory (ChRMT) j^, [g], Q, y, [q] at any fixed topological 
charge. Moreover, from the symmetry of the Dirac operator, the low-lying QCD Dirac spec- 
trum is expected to be in the same universality class of ChRMT. ChRMT thus provides a 
direct connection between Dirac eigenvalues and the effective theory describing the dynam- 
ical chiral symmetry breaking. One of the most convenient predictions of ChRMT is the 
distribution of individual eigenvalue, which can be directly compared with the lattice data. 
Such comparison has been done mainly in the quenched approximation 1^, ll,^, [l^, ex- 
cept for a work using the reweighting technique [l^ or for some recent attempts of carrying 
out dynamical fermion simulation on coarse lattices [l^. The eigenvalue spectrum in 
those calculations shows a good agreement with the prediction of ChRMT as far as the 
lattice volume is large enough ^ (1.5 fm)"^. 

In this work we perform lattice QCD simulations in and out of the e-regime including 
two light flavors of dynamical quarks. Since we are interested in the consequences of chiral 
symmetry breaking, we employ the Neuberger's overlap-Dirac operator fl?!. which pre- 



serves exact chiral symmetry 



19j at finite lattice spacings. The exact chiral symmetry is also 



helpful for numerical simulations in the e-regime, because the lowest-lying eigenvalue of the 
Hermitian overlap-Dirac operator is bounded from below (by a small but finite mass term) 
and no numerical instability occurs. The space-time volume of our lattice is x T = 16'^ x 32 
with the lattice spacing a ~ 0.11-0.125 fm. The gauge field topology is fixed to the trivial 



topological sector by introducing the extra Wilson fermions and ghosts [20]. We perform the 
Hybrid Monte Carlo simulation with the sea quark mass around 3 MeV, which corresponds 
to the e-regime: the expected pion Compton wavelength is comparable to the lattice extent 
m^L 1. The numerical cost for such a small sea quark mass is very expensive in general, 
but it is not prohibitive on the small lattice as required in the e-regime simulation. We also 
carry out simulations at several quark masses roughly in the region ms/Q-nis with the 
physical strange quark mass, which are out of the e regime. 

We study the eigenvalue spectrum of the overlap-Dirac operator on the configurations 
generated with these dynamical quarks. A good agreement of the low-lying eigenvalue 



3 



spectrum with ChRMT predictions has aheady been reported in our earher paper [2l|] for 
the run in the e-regime. The present paper describes our analysis in more detail. Since 
ChRMT provides the distribution of individual eigenvalues, the test of the agreement can 
be made using the information on the shape of the distribution, not just using the average 
values. We find a good agreement of the lowest-lying eigenvalue distribution by analyzing 
its several moments. If we look at higher eigenvalues, the agreement becomes marginal, 
because there are contaminations from the bulk of the eigenvalue spectrum corresponding 
to finite momentum pion states and other higher excited states, which are not described by 
ChRMT. We study the bulk eigenvalue spectrum and identify the region where the analysis 
in the e-regime is applied. 

A direct output from the comparison of the eigenvalue spectrum is the value of chiral 
condensate S. We extract S from the lowest-lying eigenvalue in the e-regime. For comparison 
we also calculate it on heavier quark mass lattices and extrapolate them to the chiral limit. 
Although the leading order relations in the e expansion is not valid for these lattices, the 
result in the chiral limit shows remarkable agreement with the direct calculation in the e- 
regime. We convert the value of S obtained on the lattice to the common definition in the 
continuum renormalization scheme MS using the non-perturbative renormalization (NPR) 
technique through the RI/MOM scheme which ^ is a regularization independent scheme based 
on the Green's functions of the offshell quark [22]. 

This paper is organized as follows. In Section [III ^6 review ChRMT calculations of 
the Dirac eigenvalue spectrum. The details of the numerical simulations are described in 
Section UTTl and the results of the low-lying modes in the e-regime is discussed in Section HVl 
The low-mode spectrum in the p-regime are presented in Section |Vl In Section |Vl] we also 
study the higher eigenvalue spectrum. Our conclusions are given in Section I VI II 

II. CHIRAL RANDOM MATRIX THEORY 

In the e-regime the low-lying eigenvalue spectrum of A^j-flavor QCD Dirac operator 



matches with that of Chiral Random Matrix Theory (ChRMT) 



.a 



up to a scale 



factor as described below. This can be derived by identifying the partition function of 
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V 



Nf 



(2) 



ChRMT 

^ m W 

-W^ rh 

with the QCD partition function in the e-regime. Since the dependence on the global 
topology becomes manifest in the e-regime, we work in a fixed topological sector Q. Here, 
is a complex (n + Q) x matrix, and N = 2n + Q. The parameter rh plays a role of 
quark mass. In the limit of large N, the partition function ([2]) can be modified to the form 
describing the zero-momentum mode of ChPT 



DU{det U)^ exp [— tr(mf/ + rhU^) + O(m^) 



(3) 



from which one can identify Nrh = m'EV. 



The advantage o: 
analytically known 



ChRMT ([2]) is that the eigenvalue distribution of the matrix W'^W is 
3|. Here we reproduce the known result for the case of two degenerate 
flavors and zero topological charge, which is relevant in this work. 

Let us consider the k-th. lowest microscopic eigenvalue Cfc = Nxk, with Xk the k-th. 



eigenvalue of VWHV. The distribution of (k is written as 



Cfc fik 
Jci 



Cfc- 



dCk-ii^kiCi, ■ ■ ■ ; Cfc; /^), 



(4) 



where /i = Nm = mHV. The form of uJk{Ci, ■ ■ ■ ,Ck] f^) is analytically known in the micro 
scopic limit, i.e. n oo while ^ is kept fixed: 



^kiCl, ■ ■ ■ ,Ck',fJ') 



const, e ''k/^{ 



no 



2=1 



nfc'(c? - c]y nt!(ci + det[A] ■ 



(5) 



The matrices A and B are given by 



A 



loin) /i ^hifi) 

fihifi) /o(/i) 



Bij 



CrV3(0) (3<^<A; + l) 



(1<J<2A;), (6) 



CrV4(0) {k + 2<i<2k) 



where Ci = ^Ck~ C? ^"^d fi = yCl + /^^- lii^Y^ are the modified Bessel functions. 
The spectral density is given by a sum of the individual distributions 



k 



(7) 
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FIG. 1: Low-lying spectral density in the massless limit /OrmtCCjO) (solid curve) and its decom- 
position to individual eigenvalues PkiCk',^) (dashed curves, for A; = 1, 2 and 3). The dotted curve 
represents the distribution in the infinite sea quark mass limit /9RMT(Cjoo), which corresponds to 
the quenched theory. 

In the massless and the infinite mass (or quenched) limit, it can be written in a simple form, 

PaMT(C;0) = ^ (^2^(0 - ^3(C)^i(C)) , 

PaMT(C;00) = ^(Jo'(C) + ^i'(C)), (8) 

where </j(C) denotes the Bessel functions of the first kind. Their shape and the individual 
eigenvalue distributions are shown in Figure [H 

In order to quantify the shape of the distributions, we consider n-th moments 

{Q) = JdCkC^Pk{Ck;f^), (9) 

which can be calculated numerically. The results for {{(k — (Cfe))") are shown in Figure [2] as 
a function of fi. From the plot for {(k) one can see that the lowest eigenvalue is lifted near 
the massless limit due to a repulsive force by the dynamical fermions. When /i is greater 
than 10, the eigenvalues qualitatively behave as in the quenched theory (or /i — > oo limit). 
Transition from the massless two-flavor theory to the quenched theory occurs around /i = 
1-10, where the moments of the lowest-lying eigenvalue show rather peculiar dependence on 
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FIG. 2: First (top), second (middle) and third (bottom) moments of the lowest-lying eigenvalues 
(A; = 1, 2, 3 and 4). Dependence on fi = mT,V is shown. 
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The ChRMT spectrum is expected to match with those of the QCD Dirac operator up 
to a constant T,V. For example, the lowest eigenvalue of the QCD Dirac operator Ai is 
matched as 

(Ai)/m = (Ci)/iVm = {(i)/mJ:V, (10) 

from which one can extract S, one of the fundamental constant in ChPT. Unlike the standard 
lattice QCD calculation, we do not need any chiral extrapolation, as m is already very small 
in the e-regime. By investigating the consistency with the determination through higher 
eigenvalues or their shapes, one can estimate possible systematic errors due to higher order 
effects in the e expansion. 



III. NUMERICAL SIMULATION 



A. Overlap fermion implementation 



We employ Neuberger's overlap fermion formulation [3] for the sea quarks. Its Dirac 



operator is defined as 

Dim) 



mo + 



m 



+ "^0 - 



m 



75Sgn[ifH/(-mo)], 



(11) 



where Hw = 'y^Dwi—^o) denotes the Hermitian Wilson-Dirac operator with a large neg- 
ative mass —mo- We choose mg = 1.6 throughout this work. (Here and in the following 
the parameters are given in the lattice unit.) The overlap-Dirac operator (fTTj) satisfies the 



Ginsparg- Wilson relation 



23| 



1^(0)75 + 75^(0) = —1^(0)75^(0), 
mo 



(12) 



when the quark mass m vanishes. Because of this relation, the fermion action built up with 
( ITTl) has an exact chiral symmetry under the modified chiral transformation loj ]. 

In the practical application of the overlap-Dirac operator ( ITTl) . the profile of near- zero 
modes of the kernel operator Hw{—mQ) is important, as they determine the numerical cost 
of the overlap fermion. The presence of such near-zero modes is also a problem for the 



locality property of the overlap operator 



simulations, it is known that the spectral 



24]. For most gauge actions used in practical 



density pwi^w) of the operator Hw^—mQ) is non- 
zero at vanishing eigenvalue Xw = [251] due to the so-called dislocations, i.e. local lumps 
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of the gauge configuration [26]. We avoid this problem by introducing additional fermions 
and ghosts to generate a weight 

det[Hwi-mof] 



det[Hw{—irLoy + m 



21 ' 



(13) 



in the partition function 20| . (The same idea is proposed in the context of the domain- wall 



fermion 



281].) They are unphysical as their mass is of order of lattice cutoff, and thus 
does not affect low-energy physics. The numerator suppresses the near-zero modes, while 
the denominator cancels unwanted effects for higher modes. The "twisted-mass" parameter 
rrit determines the value of threshold below which the eigenmodes are suppressed. We set 
nit = 0.2 in this work. With these extra degrees of freedom, the spectral density pwi^w) 
vanishes at the vanishing eigenvalue Xw, and the numerical cost of approximating the sign 



function in (fTTj) is substantially reduced [20 1. 



We approximate the sign function using a rational function of the form (see, e.g., |29l.l30l|) 



1 do ,,2 . ^^ bi 



iH^r ^min ^ ^ ^ + C2I-I ^ ^ 

where Xmin is the lower limit of the range of approximation and hw = Hw/Xmin- The 
coefficients 6/, q and do can be determined analytically (the Zolotarev approximation) so as 
to optimize the accuracy of the approximation. Since we have to fix the lower limit Xmin, 
we calculate a few lowest-lying eigenvalues and project them out before applying (fT4|) when 
their absolute value is smaller than Xmin- The value of Xmin is 0.144 in our simulations. The 
accuracy of the approximation improves exponentially as the number of poles n increases. 
With n = 10, the sign function sgn[Hw {—itlq)] is approximated to a 10~'^-10~^ level. Since 
the multi-shift conjugate gradient method can be used to invert all the {h^ + C2i-i)~^ terms 
at once, the numerical cost depends on n only weakly. 

In the e-regime the partition function and other physical quantities show striking depen- 
dence on the global topological charge of gauge field. With the lattice action including (fT3|) 
the topological charge never changes during the Hybrid Monte Carlo (HMC) simulations, 
which consists of molecular dynamics (MD) evolution of gauge field configuration. This is 
because the topology change must accompany a zero crossing of the eigenvalue of Hwi—fno), 
which is forbidden by the factor (|T3l) . The gauge configuration in a fixed topological sector 
can therefore be effectively sampled. In this work the simulations are restricted in the triv- 
ial topological sector Q = except for one quark mass parameter for which we carry out 



m 


traj. 


Q 


a [fm] 


0.015 


10,000 





0.1194(15) 


0.025 


10,000 





0.1206(18) 


0.035 


10,000 





0.1215(15) 


0.050 


10,000 





0.1236(14) 


0.050 


5,000 


-2 




0.050 


5,000 


-4 




0.070 


10,000 





0.1251(13) 


0.100 


10,000 





0.1272(12) 



TABLE I: Simulation parameters at (3 = 2.30. 

independent simulations at Q = —2 and —4. 

Here, we assume that the ergodicity of the simulation in a fixed topological sector is 
satisfied even with the determinant f|T3l) . In order to confirm this, we are studying the 
fiuctuation of the local topological charge density, which will be reported in a separate 
paper. 



B. HMC simulations 



We perform two-fiavor QCD simulations using the overlap fermion for the sea quarks, 
with the approximated sign function (fT^ with n = 10. Lattice size is 16^ x 32 throughout 



this work. For the gauge part of the action, we use the Iwasaki action [31|, |32| at (3 = 2.30 
and 2.35, which correspond to the lattice spacing a = 0.12 fm and 0.11 fm, respectively, 
when used with the extra Wilson fermions and ghosts. The simulation parameters are listed 
in Tables HI and HI] for (3 = 2.30 and 2.35, respectively. 

The configurations from the runs at j3 = 2.30 are for various physics measurements 
including hadron spectrum, decay constants, form factors, bag parameters, and so on. In 
this work we use them to analyze the eigenvalue spectrum. The simulation details will be 



described in a separate paper 



33| . but we reproduce some basic parameters in Tabled They 



include the sea quark mass m, trajectory length (the unit trajectory length is 0.5 MD time), 
topological charge Q and lattice spacing a determined from the Sommer scale ro (= 0.49 fm) 
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m 


traj. 


m' 


SpF2 


SpFl/SpF2 




{AH) 


p 

^ acc 


{P) 


a [fm] 


0.002 


3,690 


0.2 


0.0714 


1/4 


1/5 


0.90(23) 


0.756 


0.62482(1) 


0.1111(24) 




1,010 


0.2 


0.0625 


1/4 


1/5 


1.24(50) 


0.796 


0.62479(2) 




0.020 


1,200 


0.2 


0.0714 


1/4 


1/5 


0.035(09) 


0.902 


0.62480(1) 


0.1074(30) 


0.030 


1,200 


0.4 


0.0714 


1/4 


1/5 


0.253(20) 


0.743 


0.62480(2) 


0.1127(23) 


0.045 


1,200 


0.4 


0.0833 


1/5 


1/6 


0.189(18) 


0.768 


0.62476(2) 


0.1139(29) 


0.065 


1,200 


0.4 


0.1 


1/5 


1/6 


0.098(12) 


0.838 


0.62474(2) 


0.1175(26) 


0.090 


1,200 


0.4 


0.1 


1/5 


1/6 


0.074(19) 


0.855 


0.62472(2) 


0.1161(24) 


0.110 


1,200 


0.4 


0.1 


1/5 


1/6 


0.052(10) 


0.868 


0.62471(2) 


0.1182(22) 



TABLE II: Simulation parameters at /? = 2.35. 



34| of the heavy quark potential. In the massless hmit, the lattice spacing is found to be 
0.1184(12) fm by a linear extrapolation in m. The sea quark mass aX j3 = 2.30 covers the 
region from m^/G to with nis the physical strange quark mass. 

The runs aX (3 = 2.35 were originally intended for a basic parameter search and therefore 
the trajectory length for each sea quark mass is limited (1,200 HMC trajectories). It is 
at this (3 value that we performed a run in the e-regime by pushing the sea quark mass 
very close to the chiral limit m = 0.002, which is one order of magnitude smaller than the 
sea quark mass in other runs. In Table [TTl we summarize several simulation parameters. 
Among them, the basic parameters are the sea quark mass m, trajectory length, plaquette 
expectation value (P), and lattice spacing. The massless limit of the lattice spacing is 
evaluated to be 0.1091(23) fm using a linear extrapolation with data above m = 0.020. This 
value is consistent with the result of the e-regime run at m = 0.002. The other parameters 
are explained below. 

The HMC simulation with the overlap fermion was first attempted by Fodor, Katz and 
Szabo [35] and soon followed by two other groups They introduced the so-called 

reflection-refraction trick in order to treat the discreteness of the HMC Hamiltonian at 
the topological boundary. This leads to a significant additional cost for dynamical overlap 
fermions compared to other (chirally non-symmetric) fermion formulations. We avoid such 
extra costs by introducing the extra Wilson fermion determinants ( |T3l) . with which the MD 
evolution never reaches the topological boundary. 



11 



In the implementation of the HMC algorithm, we introduce the Hasenbusch's mass pre- 
conditioner 38|] together with the multiple time step technique 39|]. Namely, we rewrite the 
fermion determinant as 



det[D{m)Y = det[D(m')]^det 



Dim) 



Dim 



>\2 



(15) 



by introducing a heavier overlap fermion with mass m' . We then introduce a pseudo-fermion 
field for each determinant. In the right hand side of ( |T5l) the second term is most costly as 
it requires an inversion of the overlap operator with a small mass m. On the other hand, 
the contribution to the MD force from that term can be made small by tuning m' close 
to m. With the multiple time step technique, such small contribution does not have to be 
calculated frequently, while the force from the first term must be calculated more often. 
We introduce three time steps: (i) 5tpp2 for the ratio det[D(m)^/D(m')^], (ii) drppi for 
the preconditioner det[i5(m')]^, and (iii) 5tg for the gauge action and the extra Wilson 
fermions ( |T3i) . By investigating the size of MD forces from each term, we determine the 
time steps and the preconditioner mass m' as listed in Table [Tll For the run in the e- 
regime {(5 = 2.35, m = 0.002) we switched Stpp2 to a smaller value in the middle of the 
run, since we encounter a trajectory which has exceptionally large MD force from the ratio 
det[D{my / D{m'y] probably due to a small eigenvalue of D{m). 

An average shift of Hamiltonian during a unit trajectory (AH) determines the acceptance 
rate Pace in the HMC algorithm. It must be 0(1) or less to achieve a good acceptance rate, 
which is satisfied in our runs as listed in Table [Tll The value at m = 0.002 is larger and 
around 0.9-1.2. This is due to so-called "spikes" phenomena, i.e. exceptionally large values 
(~ 0(10 — 100)) of AH at some trajectories. The spikes are potentially dangerous as they 
may spoil the exactness of the HMC algorithm, but we believe that this particular run is 
valid since we have checked that the area preserving condition {e~^^) = 1 is satisfied within 
statistical errors. 

For the inversion of the overlap operator we use the relaxed conjugate gradient algorithm 
[4^. The trick is to relax the convergence condition of the inner solver as the conjugate 
gradient loop proceeds. This is allowed because the change of the solution vector becomes 
smaller at the later stages of the conjugate gradient. The gain is about a factor of 2 compared 
to the conventional conjugate gradient. In the middle of the simulations at /3 = 2.30, we 



replaced the overlap solver by the one with a five-dimensional implementation 
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4l|. This is 
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FIG. 3: Number of the Wilson-Dirac operator multiplication per trajectory (upper panel) and per 
an overlap inversion (lower panel) for (5 = 2.35. The curves are fits to data above m = 0.030 with 
the form oc 1/m". 

faster by another factor of 4-5 than the relaxed conjugate gradient method. These details 
of the algorithm will be discussed in a separate paper [SSj. 

The numerical cost depends on how precisely the matrix inversions are calculated. At 
an inner level there are inversions of the Hermitian Wilson-Dirac operator appearing in the 
rational approximation ^\M . The n inversions can be done at the same time using the multi- 
shift conjugate gradient. We calculate until all the solutions reach the relative precision 10~^ 
when adopted in the calculation of the HMC Hamiltonian. This value matches the precision 
we are aiming at for the approximation of the sign function. In the molecular dynamics 
steps the relative precision is relaxed to 10^^. The conjugate gradient for the overlap-Dirac 
operator at the outer level is also carried out to the level of the 10^^ (10^'') relative precision 
in the HMC Hamiltonian (MD force) calculation. 

The numerical cost can be measured by counting the number of the Wilson-Dirac operator 
multiplication, although other manipulations, such as the linear algebra of vectors, are not 
negligible. The number of the Wilson-Dirac operator multiplication is plotted in Figure [3] 
for the runs aX (3 = 2.35. The upper panel shows the cost per trajectory; the lower panel 
presents the cost of inverting the overlap-Dirac operator when we calculate the Hamiltonian 
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at the end of each trajectory. The expected mass dependence for the overlap solver is 
1/ \Jm? + |AiP with Ai the lowest-lying eigenvalue of the overlap operator -D(O). Therefore, 
the cost is proportional to 1/m only when m is much greater than |Ai|. This condition is 
satisfied for m at and larger than 0.030, where |Ai| is around 0.004 as we show later. Fitting 
the data with the scaling law ~ l/m° above m = 0.030, we obtain the power a as 0.82, 
which is roughly consistent with the expectation. For the total cost of the HMC Hamiltonian 
(upper panel), the quark mass dependence is more significant, since it depends on the choice 
of the step sizes. It is not even a smooth function of m. If we fit the data with the power 
law ~ above m = 0.030 as in the case of the solver, we obtain a = 0.49, which gives 

a much milder quark mass dependence. 

The machine time we spent is roughly one hour per trajectory for the run in the e-regime 
(m = 0.002) on a half rack (512 computing nodes) of IBM BlueGene/L. The cost at other 
mass parameters is lower as one can see in Figure [31 The numerical cost at /5 = 2.30 is 
higher, because the number of the near-zero modes of ifvy(— mo) is significantly larger. 

For comparison we also generated quenched configurations on a 16^ x 32 lattice at j3 = 
2.37 in the topological sector Q = and 2. We must use the HMC algorithm even for the 
quenched simulation, as it contains the extra Wilson fermions f fT3|) . We accumulated 20,000 
trajectories for each topological sector and used the gauge configurations for measurement 
at every 200 trajectories. The lattice spacing is 0.126(2) fm, which matches the dynamical 
lattices at /5 = 2.30 in the heavier sea quark mass region m = 0.075 and 0.100. In the chiral 
limit the dynamical lattices are slightly finer. 

C. Eigenvalue calculation 

In the HMC simulations described in the previous section, we stored the gauge config- 
urations at every 10 trajectories for measurements. For those configurations we calculate 
lowest 50 eigenvalues and eigenvectors of the overlap-Dirac operator -D(O). In the analysis 
of this work we only use the eigenvalues. 

We use the implicitly restarted Lanczos algorithm for a chirally projected operator 

D+ = P+I)(0)P+, (16) 
where P+ = (1 + 75)/2. This operator is Hermitian and its eigenvalue gives the real part of 
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FIG. 4: Ensemble averages of the lowest five eigenvalues {Xk) {k = 1-5) as a function of sea quark 
mass at /? = 2.35. Dashed line shows A = m. 

the eigenvalue of the original overlap operator D{0). The pair of eigenvalues A°'' (and its 
complex conjugate) of -D(O) can be obtained from ReA°" using the relation |1 — A°^/mop = 1 
derived from the Ginsparg- Wilson relation (|T2l) . 

In the calculation of the eigenvalues we enforce better accuracy in the approximation 
of the sign function by increasing the number of poles in the rational function. The sign 
function is then approximated at least to the 10~^^ level. In order to improve the convergence 
of the Lanczos algorithm we use the Chebyshev acceleration technique 42, 4^ and optimize 
the window of eigenvalues for the target low-lying modes. 

For the comparison with ChRMT, the lattice eigenvalue X""" is projected onto the imagi- 
nary axis as A = ImA°^/(l — ReA°^/(2mo)). Note that A is very close to ImA°^ (within 0.05%) 
for the low-lying modes we are interested in. We consider positive A's in the following. 

In Figure H] we plot the ensemble averages of the lowest 5 eigenvalues (A^) {k = 1-5) as 
a function of the sea quark mass. The data at (3 = 2.35 are shown. We observe that the 
low-lying spectrum is lifted as the chiral limit is approached. This is a direct consequence 
of the fermion determinant ~ rifod-^feP + fn"^), which repels the small eigenvalues from zero 
when the lowest eigenvalue is larger than m. This is exactly the region where the numerical 
cost saturates as it is controlled by Ai rather than m. 

Figure [5] shows a Monte Carlo history of the lowest-lying eigenvalue Ai at the lightest (m 
= 0.002) and the heaviest (m = 0.110) sea quark masses at /5 = 2.35. At m = 0.002 we find 
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FIG. 5: Monte Carlo history of the lowest eigenvalue Ai for the sea quark masses m = 0.002 (top) 
and 0.110 (bottom) at /3 = 2.35. 
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FIG. 6: Jackknife bin-size dependence of the error for the eigenvalue average (A^) (k = 1-4) at 
P = 2.35 and m = 0.002. 

some long range correlation extending over a few hundred trajectories, while the history m 
= 0.110 seems more random. In order to quantify the effect of autocorrelation we investigate 
the bin-size dependence of the jackknife error for the average (A^) {k = 1-5). As can be seen 
from Figure [6]the jackknife error saturates around the bin-size 20, which corresponds to 200 
HMC trajectories. This coincides with our rough estimate from Figure O In the following 
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analysis we take the bin-size to be 20 at m = 0.002 and 10 at other sea quark masses. 



IV. LOW-MODE SPECTRUM IN THE e-REGIME 

In this section we describe a comparison of the lattice data for the low-lying eigenvalues 
with the predictions of ChRMT. The most relevant data set in our simulations is the one at 
m = 0.002 and (3 = 2.35, since this is the only run within the e-regime. 

First we determine the scale, or the chiral condensate, from the first eigenvalue through 
f lTU]) . By solving 



recursively in order to correct the fi dependence of (Ci), we obtain /i = 0.556(16) and 
j^iat _ 0.00212(6) in the lattice unit. In the physical unit, the result corresponds to S'*^* = 
[240(2) (6) MeV]^ where the second error comes from the uncertainty in the lattice scale a 
= 0.107(3) fm. In the above, we put a superscript '/at' to the chiral condensate S in order 
to emphasize that it is defined on the lattice. The error of {(i) = 4.30 from the statistical 
error of (Ai) is neglected (within 0.1%). Note that /i = 0.556 is already very close to the 
chiral limit as one can see from Figure O For the average of the lowest eigenvalue (Ci) the 
difference from the massless limit is only 0.9%. 

Next, let us compare the higher eigenvalues of the Dirac operator. We plot the ratios 
(Cfc) / (0) of eigenvalues in Figure [71 The lattice data agree well with the ChRMT predictions 
(middle panel). It is known that there exists the so-called flavor-topology duality in ChRMT: 
the low-mode spectrum is identical between the two-flavor (massless) theory at Q = and 
the quenched theory at \Q\ = 2 (right panel), while the quenched spectrum at Q = is 
drastically different (left panel). This is nicely reproduced by the lattice data. Note that 
the finite /i(~ 0.56) corrections to the massless case are very small. 

Another non-trivial comparison can be made through the shape of the eigenvalue distri- 
butions. We plot the cumulative distribution 



of the three lowest eigenvalues in Figure El The agreement between the lattice data and 
ChRMT (solid curves) is quite good for the lowest eigenvalue, while for the higher modes the 
agreement is marginal. This observation can be made more quantitative by analyzing the 
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(Ci)//^, fi = mEV, 
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FIG. 7: Ratio of the eigenvalues {Ck)/{Ci) ^or combinations of k and / G 1~4 (denoted in the plot 
as k/l). We use the input, /i = 0.556(16), which is obtained from the lowest eigenvalue average. 
In addition to the two- flavor QCD data (middle), quenched data at \Q\ =0 (left) and 2 (right) 
at P = 2.37 are shown. Lattice data (circles) are compared with the ChRMT predictions (bars). 
Note that the finite /i(~ 0.56) corrections to the massless case are tiny. 
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FIG. 8: The accumulated histogram of the eigenvalues, x-error comes from the statistical error 
of S. The solid lines are the ChRMT results with an input for T, from the average of the lowest 
eigenvalue. 
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1 4.30 [4.30] 


1.52 


1.48(12) 


0.41 


0.74(27) 


2 7.62 7.25(13) 


1.73 


2.11(24) 


0.28 


0.83(43) 


3 10.83 9.88(21) 


1.88 


2.52(31) 


0.22 


0.38(58) 


4 14.01 12.58(28) 


2.00 


2.39(31) 


0.18 


0.22(66) 



TABLE III: Moments of the low-lying eigenvalues. Comparison between CliRMT and lattice data 
are made for the first three moments. The average value of the lowest eigenvalue (Ci) = {Xi)T,V is 
an input for S. Here, the errors of (Cfe)'s or their higher moments due to the uncertainty of S are 
neglected (within 0.1%). 

moments defined in (Q. In Table UTTl we list the numerical results of both ChRMT and lattice 
data for the subtracted moments {{Ck — (Cfc))")- The overall agreement is remarkable, though 
we see deviations of about 10% in the averages. The deviations in the higher moments are 
larger in magnitude but statistically less significant (less than two standard deviations). 

The leading systematic error in the determination of S is the finite size effect, which scales 
as O(e^) ~ 0{1/{Ft,LY). Unfortunately we can not calculate such a higher order effect 
within the framework of ChRMT, but we can estimate the size of the possible correction 
using the higher order calculations of related quantities in ChPT. To the one-loop order, the 
chiral condensate is written as 

N}-1 A 



1 + 



Nf (F^L) 



(19) 



where jSi is a numerical constant depending on the lattice geometry [44,]- The value for the 
case of the x (2L) lattice is 0.0836. Numerically, the correction is 13% assuming the pion 
decay constant to be = 93 MeV. 

The most direct way of reducing the systematic error is to increase the volume, which is 
very costly, though. Other possibility is to check the results with quantities for which the 
higher order corrections are known. Meson two-point functions in the e-regime are examples 
of such quantities. A work is in progress to calculate the two-point functions on our gauge 
ensembles. 

We quote the result of E in the continuum regularization scheme, i. e. the MS scheme. We 
have calculated the renormalization factor Zg^^{2 GeV) using the non-perturbative renor- 



malization technique through the RI/MOM scheme 
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22 1 . Calculation is done on the e- 



regime (m = 0.002) lattice with several different valence quark masses. The result is 
Zg^^{2 GeV) = 1.14(2). Details of this calculation will be presented in a separate paper. 
Including the renormalization factor, our result is 

E^(2 GeV) = [251(7)(11) MeV]l (20) 

The errors represent a combined statistical error (from Ai, tq, and Z^^{2 GeV)) and the 
systematic error estimated from the higher order effects in the e-expansion as discussed 
above. Since the calculation is done at a single lattice spacing, the discretization error 
cannot be quantified reliably, but we do not expect much larger error because our lattice 
action is free from 0(a) discretization effects. 

V. LOW-MODE SPECTRUM IN THE j5-REGIME 

For heavier sea quarks, the e-expansion is not justified and the conventional p-expansion 
should be applied instead. Therefore, the correspondence between the Dirac eigenvalue 
spectrum and ChRMT is not obvious. On the other hand, for heavy enough sea quarks the 
low-lying eigenvalues should behave as if they are in the quenched lattices. Here we assume 
that the correspondence is valid in the intermediate sea quark mass region too, and compare 
the lattice data with the ChRMT predictions for larger /i = rnSV . Strictly speaking, the 
theoretical connection to ChRMT is established only at the leading order of the e expansion, 
which is vahd when (M^L)^ ~ (mSV^)/(F^L)^ <^ 1 is satisfied. 

In Figure [9] we plot the eigenvalue ratios (Afc) / (Ai) {k = 2-4) as a function of m'EV. The 
data are shown for both f3 = 2.35 and 2.30. The curves in the plots show the predictions 
of ChRMT. The expected transition from the dynamical to quenched lattices can be seen 
in the lattice data below mTiV ~ 10. The mass dependence at (3 = 2.35 is consistent with 
ChRMT within relatively large statistical errors, while the precise data at (3 = 2.30 show 
some disagreement especially for third and fourth eigenvalues. 

We extract the chiral condensate S for each sea quark mass using the same method 
applied in the e-regime taking account of the mass dependence of {(i). The results at f3 = 
2.30 are plotted in Figure [10] (open circles). We use a physical unit for both m and S''^*; the 
lattice scale is determined through ro after extrapolating the chiral limit. The results show 
a significant sea quark mass dependence. If we extrapolate linearly in sea quark mass using 
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FIG. 9: Sea quark mass dependence of the ratio of the eigenvalues (Afe)/(Ai) for k = 2, 3, and 4. 
Data at /? = 2.35 (top) and 2.30 (bottom) are shown. Horizontal error comes from the uncertainties 
of T, obtained in the e-regime. The quenched results at /3 = 2.37 with Q = (left) and Q = 2 
(right) are also plotted to see the flavor-topology duality. 

three lowest data points we obtain S'"* = [245(5) (6) MeV]^ in the chiral limit. This value is 
consistent with the result in the e-regime as shown in the plot. 

In Figure dn] we also plot data points for non-zero topological charge {\Q\ = 2 and 4) at 
m = 0.050. We find some discrepancy between \Q\ = and 2 while \Q\ = 4 is consistent 
with \Q\ = 0. The size of the disagreement is about 4% for (s'"*)V3 g^^^ ^^j^^g 12% for E'"*, 
which is consistent with our estimate of the higher order effect in the e expansion. 
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FIG. 10: Sea quark mass dependence of the chiral condensate (S'"*)^/^ extracted from the lowest 
eigenvalue. Open symbols denote the data at (3 = 2.30 with their chiral extrapolation shown by a 
filled circle. A filled square is the result in the e-regime (/3 = 2.35 and ma = 0.002). The lattice 
scale is determined through the chiral extrapolation of rg; its statistical error is not taken into 
account in the plot. 

VI. BULK SPECTRUM 



Although our data for the Dirac eigenvalue spectrum show a qualitative agreement with 
the ChRMT predictions, there are 0(10%) deviations, which is significant for the larger 
eigenvalues as seen in Table. IIIII This can be understood by looking at higher eigenvalue 
histogram, which we call the bulk spectrum. Figure [TT] shows a histogram of 50 lowest 
eigenvalues in the e-regime {j3 = 2.35, m = 0.002). The normalization is fixed such that it 
corresponds to the spectral density 



k 



(21) 



divided by the volume in the limit of vanishing bin size. 

In order to understand the shape of the data in Figure [TT] at least qualitatively, we 
consider a simple model. Away from the low-mode region one expects a growth of the 
spectral function as ~ 3A'^/47r^, which is obtained from the number of plain- wave modes of 
quarks in the free case. By adding the condensate contribution S/tt from the Banks-Casher 
relation 45|] we plot a dashed curve in Figure [Til Near the microscopic limit XT,V 0, 
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FIG. 11: Eigenvalue histogram of the lowest 50 eigenmodes. The bold curve shows the ChRMT 
prediction of the spectral density and the dashed line is (free theory + constant S/vr), in which we 
use S = 0.00212 obtained in the e-regime. 

the ChRMT prediction Sprmt(AS1^; rriT^V) is expected to match with the data, where Prmt 
is defined by (I7j). We plot the massless case T,p^mt{^^V]0) in Figure fTTl for a comparison. 
(Deviation of the spectrum at mSy = 0.56 from the massless case is only ~ 1%.). 

The ChRMT curve gives a detailed description of the Banks-Casher relation: it ap- 
proaches a constant E/vr in the large volume limit. On the other hand, since ChRMT is 
valid only at the leading order of the e expansion, the region of 0{X^) growth cannot be 
described. Therefore, for the analysis of the microscopic eigenvalues to be reliable, one has 
to work in a flat region where the O(A^) contribution is negligible. This is the reason that 
the lowest eigenvalue is most reliable to extract S in our analysis in the previous sections. 

From Figure [TT] we observe that the flat region does not extend over ASV ~ 15, which 
roughly corresponds to the fourth lowest eigenvalue in our data. Already at around this 
upper limit, the eigenvalues are pushed from above by a repulsive force from the bulk eigen- 
modes rapidly increasing as oc A^, and the ratio (Afc)/(Ai) is systematically underestimated 
for k = 3 and 4 as found in Figure [71 This effect is regarded as one of the finite size effect, 
because the A^ term scales as (ASV^)^/(SV^)^ and its magnitude in the microscopic regime 
is suppressed for larger volumes as In addition, the peaks of the first few eigenvalues 

move towards ASV = for larger volumes, and thus become less sensitive to the effects from 
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FIG. 12: Eigenvalue histogram for the (5 = 2.30 lattices. Solid curves show quenched ChRMT 
and asymptotic form obtained from the free quark theory. For a normalization we use S = 0.00212 
obtained in the e-regime. 

bulk eigenmodes. 

The bulk spectrum for heavier sea quark masses, which are out of the e-regime, is also 
interesting in order to see what happens after the transition to the "quenched-like" region 
of the eigenvalue spectrum. In Figure [12] the eigenvalue histogram is shown for (3 = 2.30 
lattices at m = 0.015, 0.035, 0.050 and 0.070, all of which are in the p-regime. The plot 
is normalized with S = 0.00212, which is the value after the chiral extrapolation shown 
in Figure [10 First of all, the physical volume at /3 = 2.30 is about 30% larger than that 
at /? = 2.35. Therefore, as explained above, the growth of 0{\^) is expected to be much 
milder and the lattice data is consistent with this picture. The flat region extends up to 
around mSy ~ 30. Second, because the microscopic eigenvalue distribution approaches 
that of the quenched theory, the lowest peak is shifted towards the left. Overall, the number 
of eigenvalues in the microscopic region increases a lot. Unfortunately, the correspondence 
between ChPT and ChRMT is theoretically less clear, since the sea quark masses are in 
the p-regime. In order to describe this region, the standard ChPT must be extended to the 
partially quenched ChPT and a mixed expansion has to be considered. Namely, the sea 
quarks are treated in the p-expansion, while the valence quarks are put in the e-regime to 
allow the link to ChRMT. In this paper we simply assume that ChRMT can be applied for 
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finite sea quark masses out of tlie e-regime. We observe in Figure [12] tliat the distribution 
near tlie lowest eigenvalue is well described by CliRMT, but the peak grows as the quark 
mass increases. This means that the effective value of S grows as the quark mass increases, 
which is consistent with the sea quark mass dependence of S plotted in Figure (TDl 

VII. CONCLUSIONS 

We studied the eigenvalue spectrum of the overlap-Dirac operator on the lattices with 
two-flavors of dynamical quarks. We performed dynamical fermion simulation in the e- 
regime by pushing the sea quark mass down to 3 MeV. For comparison, we also calculated 
the eigenvalue spectrum on the p-regime lattices at two lattice spacings with sea quark mass 
in the range m^/G-ms. All the runs are confined in a fixed topological charge Q = 0, except 
for a few cases with finite Q. 

We found a good agreement of the distribution of low-lying eigenvalues in the e-regime 
with the predictions of ChRMT, which implies a strong evidence of the spontaneous breaking 
of chiral symmetry in Nf = 2 QCD. We extracted the chiral condensate as S'^^(2 GeV) = 
[251(7)(11) MeV]'^ from the lowest eigenvalue. The renormalization factor was calculated 
non-perturbatively. The value of S contains a systematic error of ~ 10% due to the higher 
order effect in the e expansion 0{1/ Ft^L). Better determination of S will require larger 
physical volumes to suppress such finite size effects. 

Out of the e-regime (the case with heavier sea quark masses) the Dirac eigenvalue distri- 
bution still shows a reasonable agreement with ChRMT. The value of S extracted in this 
region shows a significant quark mass dependence, while its chiral limit is consistent with 
the e-regime result. 

Further information on the low-energy constants can be extracted in the e-regime by 
calculating two- and three-point functions or analyzing the Dirac eigenvalue spectrum with 



imaginary chemical potential The present work is a first step towards such 

programs. 
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